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We discuss the underlying relativistic physics which causes 
neutron stars to compress and collapse in close binary systems 
as has recently been observed in numerical (3+1) dimensional 
general relativistic hydrodynamic simulations. We show that 
compression is driven by velocity-dependent relativistic hy- 
drodynamic terms which increase the self gravity of the stars. 
They also produce fluid motion with respect to the corotat- 
ing frame of the binary. We present numerical and analytic 
results which confirm that such terms are insignificant for uni- 
form translation or when the hydrodynamics is constrained to 
rigid corotation. However, when the hydrodynamics is uncon- 
strained, the neutron star fluid relaxes to a compressed non- 
synchronized state of almost no net intrinsic spin with respect 
to a distant observer. We also show that tidal decompression 
effects are much less than the velocity-dependent compression 
terms. We discuss why several recent attempts to analyze 
this effect with constrained hydrodynamics or an analysis of 
tidal forces do not observe compression. We argue that an 
independent test of this must include unconstrained relativis- 
tic hydrodynamics to sufficiently high order that all relevant 
velocity-dependent terms and their possible cancellations are 
included. 
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I. INTRODUCTION 

The physical processes occurring during the last orbits 
of a neutron-star binary are currently a subject of intense 
interest [Q- Qj. In part, this recent surge in interest 
stems from relativistic numerical hydrodynamic simula- 
tions in which it has been noted |lf|| that as the stars 
approach each other their interior density increases. In- 
deed, for an appropriate equation of state, our numerical 
simulations indicate that binary neutron stars collapse 
individually toward black holes many seconds prior to 
merger. This compression effect would have a significant 
impact on the anticipated gravity- wave signal from merg- 
ing neutron stars. It could also provide an energy source 
for cosmological gamma-ray bursts . 

In view of the unexpected nature of this neutron star 
compression effect and its possible repercussions, as well 
as the extreme complexity of strong field general rela- 
tivistic hydrodynamics, it is of course imperative that 



there be an independent confirmation of the existence of 
neutron star compression before one can be convinced 
of its operation in binary systems. In view of this it is 
of concern that the initial numerical results reported in 
jl]-[| have been called into question. A number of re- 
cent papers P~p^| have not observed this effect in New- 
tonian tidal forces [Q, first post-Newtonian (IPN) dy- 
namics P~ppd[|, tidal expansions |L0-12|, or in binaries 
in which rigid corotation has been imposed []l3|| . The 
purpose of this paper is to point out that none of these 
recent studies could or should have observed the com- 
pression effect which we observe in our calculations. 

Moreover, this flurry of activity has caused some con- 
fusion as to the physics to which we attribute the effects 
observed in the numerical calculations. The present pa- 
per, therefore, summarizes our derivation of the physics 
which drives the collapse. We illustrate how such terms 
have been absent in some Newtonian or post-Newtonian 
approximations to the dynamics of the binary system. 
We also present numerical results and analytic expres- 
sions which demonstrate how the compression forces re- 
sult in an orbiting dynamical system from the presence 
of fluid motion with respect to the corotating frame. As 
such, they could not appear in an an analysis of relativis- 
tic external tidal forces no matter how many orders are 
included in the tidal expansion parameter (e.g. JiT| , p^[ ) 
unless self gravity from internal hydrodynamic motion is 
explicitly accounted for. The effect could also not arise 
in systems with uniform translation or rigid corotation. 

The implication of the present study is that any at- 
tempt to confirm or deny the compression driving force 
requires an unconstrained, untruncated relativistic hy- 
drodynamic treatment. At present, ours is still the only 
existing such calculation. Hence, despite claims to the 
contrary ||-[l4|], the neutron star compression effect has 
not yet been independently tested. 

Another confusing aspect surrounding the numerical 
results has been our choice of a conformally flat spatial 
three-metric for the solution of the field equations. In- 
deed, it has been speculated that this approximate gauge 
choice (in which the gravitational radiation is not explic- 
itly manifest) may have somehow led to spurious results. 
A second purpose of this paper, therefore, is to empha- 
size that the compression driving terms are a completely 
general result from the relativistic hydrodynamic equa- 
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tions of motion. The advantages of the conformally flat 
condition are that the algebraic form of the compression 
driving terms is easier to identify and that the solutions 
to the field equations obtain a simple form. It does not 
appear to be the case, however, that the imposition of 
a conformally flat metric drives the compression. It has 
been nicely demonstrated in the work of Baumgarte et 
al. jl3| that conformal flatness does not necessarily lead 
to neutron-star compression. 



II. THE SPATIALLY CONFORMALLY FLAT 
CONDITION 

There has been some confusion in the literature as to 
the uncertainties introduced by imposing a conformally 
flat condition (henceforth CFC) on the spatial three- 
metric. Therefore we summarize here some attempts 
which we and others have made to quantify the nature 
of this approximation. 

The only existing strong field numerical relativistic hy- 
drodynamics results in three unrestricted spatial dimen- 
sions to date have been derived in the context of the CFC 
as described in detail in Jl]-|3|]. 

We begin with the usual ADM (3+1) metric 017) 
in which there is a slicing of the spacetime into a one- 
parameter family of three-dimensional hypersurfaces sep- 
arated by differential displacements in a timelike coordi- 
nate, 



ds 2 = -(a 2 - Pifi^dt 2 + 2j3 i dx i dt + j ij dx i dx j 



(1) 



where we take Latin indices to run over spatial coordi- 
nates and Greek indices to run over four coordinates. We 
also utilize gcomctrizcd units (G = c = 1) unless other- 
wise noted. The scalar a is called the lapse function, [3i 
is the shift vector, and jij is the spatial three metric. 

In what follows, we make use of the general relation 
between the determinant of the four metric g a p and the 
ADM metric coefficients 



det(g a p) = -a 2 dei(7ij) 



2 2 

a 7 



(2) 



where 7 = ^J—det(^ij). 

The conformally flat metric condition simply expresses 
the three metric of Eq. ([[]) as a position dependent con- 
formal factor 4> 4 times a flat-space Kronecker delta 



jij = (j> 4 5i. 



(3) 



It is common practice (e.g. [|18| f20|) to impose this 
condition when solving the initial value problem in nu- 
merical relativity. It is the natural choice for our three- 
dimensional quasiequilibrium orbit calculations j^j which 
in essence seek to identify a sequence of initial data con- 
figurations for neutron-star binaries. 

The reason conformal flatness is chosen most fre- 
quently for the initial value problem is that it simplifies 
the solution of the hydrodynamics and field equations. 



The six independent components of the three metric are 
reduced to a single position dependent conformal factor. 

Since conformal flatness implies no transverse traceless 
part of 7ij it can minimize the amount of initial grav- 
itational radiation apparent in the initial configuration. 
However, in general the physical data still contain a small 
amount of preexisting gravitational radiation. This has 
been clearly demonstrated in numerical calculations of 
axisymmetric black- hole collisions pl| ]. In exact numer- 
ical simulations, the gravitational radiation appears as 
the time derivatives of the spatial three metric (7^ ) and 
its conjugate (the extrinsic curvature Kij) are evolved. 
The immediate evolution of the fields from conformally 
flat initial data is characterized by the development of a 
weak gravity wave exiting the system. 

An estimate of the radiation content of initial data 
slices for axisymmetric black hole collisions has been 
made by Abrahams [p2) . Even for high values of momen- 
tum, the initial slice radiation is always less than about 
10% of the maximum possible radiation energy (as esti- 
mated from the area theorem) . 

Two questions then are relevant to our application of 
the CFC. One is the validity of this metric choice for 
the initial value problem, and the other is the effect on 
the system of the "hidden" gravitational radiation in the 
physical data. 

Regarding the validity of the CFC one has a great deal 
of freedom in choosing coordinates and initial conditions 
as long as the initial space is Riemannian and the met- 
ric coefficients satisfy the constraint equations of general 
relativity pSf] . Indeed, we have shown in (2) that ex- 
act solutions for the CFC metric coefficients can be ob- 
tained by imposing the ADM Hamiltonian and momen- 
tum constraint conditions. Nevertheless, in three dimen- 
sions a physical space is conformally flat if and only if 
the Cotton- York tensor vanishes |25|,E!6| , 



C ij = 2ei ki 1 R j _ Kj R 



(4) 



where R J k is the Ricci tensor and R is the Ricci scalar 
for the three space. 

Equation ([|) vanishes by fiat for the three-space met- 
ric we have chosen. However, conformally flat solu- 
tions for physical problems have only been proven [^5 26 
for spaces of special symmetry (e.g. constant curvature, 
spherical symmetry, time symmetry, Robertson- Walker, 
etc. [E5B). Hence, the invocation of the CFC here and 
in other applications is an assumption. That is, it is a 
valid solution to the Einstein constraint equations, but 
does not necessarily describe a physical configuration to 
which two neutron stars will evolve. Nevertheless, this is 
a valid approximation as long as the nonconformal con- 
tributions from the 7^ and Kij equations in the exact 
two-neutron star problem remain small. Indeed, numer- 
ical tests for an axisymmetric rotating neutron star |27j] 
and a comparison of the CFC vs. an exact metric ex- 
pansion for an equal-mass binary Q have indicated that 
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conformal flatness is a good approximation when it can 
be tested. 

As a related illustration, consider the Kerr solution for 
a rotating black hole. It is well known that the Kerr met- 
ric is not conformally flat. The close binaries we study 
have specific angular momentum only slightly greater 
than that of an extreme Kerr black hole. Also, they ul- 
timately merge and collapse to a single Kerr black hole. 
Hence, an analysis of the Cotton- York tensor for a Kerr 
black hole is another indicator of the degree to which con- 
formal flatness is a valid approximation for neutron-star 
binaries. 

Figure |l|. gives the dimensionless scaled Cotton- York 
parameter C e ^m 3 for a maximally rotating Kerr black 
hole as a function of proper distance. For illustration, 
consider the decrease of this quantity as one moves away 
from the horizon at to = r as a measure of the rate at 
which the metric becomes conformally flat. The maxi- 
mally rotating (a = to) black hole of this example, how- 
ever, is an extreme example of compactness and angular 
velocity relative to any of the neutron stars in our simu- 
lations. 

It can be seen in figure [l] that, even for this extreme 
case, the dimensionless tensor coefficient C e ^rv? dimin- 
ishes rapidly away from the black hole. At the separation 
of interest for binary neutron stars approaching their fi- 
nal orbits (r/m ~ 25 where to is the total binary mass 
and r the separation between stars), this coefficient has 
already diminished to ~ I0~ 3 of the value at the event 
horizon, (r/m ~ I). Thus, the effect of either star on 
its companion is probably well approximated by confor- 
mal flatness. Regarding the interior of the neutron stars 
themselves, in our studies the stars are rotating so slowly 
(even when corotating) that the deviation from conformal 
flatness is probably negligible. Thus, it seems plausible 
that conformal flatness is a reasonable approximation for 
most physical aspects involving the spatial three-metric 
of binary neutron-star systems. 

The next issue concerns the "hidden" radiation in the 
physical data. To address this we decompose the extrin- 
sic curvature into longitudinal KV and transverse 
components as proposed by York |2q ] , 

A"' = K% + K% . (5) 

By definition the transverse part obeys 

DiK$=0 , (6) 

where Di are covariant derivatives. The longitudinal part 
can be derived from a properly symmetrized vector po- 
tential. We find 

DiK^SnS* , (7) 

where S l are spatial components of the contravariant 
four-momentum density. 

The product Kjl Kxij is a measure of the hidden radia- 
tion energy density. To find Kjl then from our numerical 



calculations, we first find Kij by choosing maximal slic- 
ing [Tr(Kij) — 0] and requiring that the trace free part 
of the jij equation vanish. This gives B 

2aK ij = (D i p j +D j /3 i -^ct>-%D k p k ) . (8) 

We then determine K l £ from the equilibrium momentum 
density [Eq. (0)] and subtract K% from K ij . 

We find that this measure of the " hidden" gravitational 
radiation energy density is a small fraction of the total 
gravitational mass energy of the system, 

/ K$ 'K Ti ~ a 2 x 1(T 5 M G . (9) 

J OTT 

Hence, we conclude that the CFC is probably a good 
approximation to the initial data. 

This should be an excellent approximation for the de- 
termination of stellar structure and stability. However, 
an unknown uncertainty enters if one attempts to recon- 
struct the time evolution of the system (e.g. the gravita- 
tional waveform) from this sequence of quasistatic initial 
conditions. At present we make this connection approx- 
imately via a multipole expansion |24f| for the gravita- 
tional radiation as described in [El. 



A. An Electromagnetic Analogy 

The meaning of imposing a conformally flat spatial 
metric can, perhaps, be qualitatively understood in an 
electromagnetic analogy. Both the ADM formulation 
of relativity and Maxwell's equations can be written as 
two constraint equations plus two dynamical equations. 
In electromagnetism the constraint equations for electric 
and magnetic fields are embodied in the V • E and V • B 
equations, while the dynamical equations are contained 
in Ampere's law and Faraday's law. In relativity the 
analogous constraint equations are the ADM momentum 
and Hamiltonian constraints. The dynamical equations 
are the ADM and 7y equations. In either electro- 
magnetism or gravity, any field configuration which satis- 
fies the constraint equations alone represents a valid ini- 
tial value solution. However, one must analyze its phys- 
ical meaning. 

For example, consider two orbiting charges. One could 
construct an electric field which satisfies the constraint 
by simply summing over the electrostatic field from two 
point charges. Similarly, one can construct a static mag- 
netic field from the charge current by imposing E = B = 
in the dynamical equations. However, by forcing the 
dynamical equations to vanish, one has precluded the ex- 
istence of electromagnetic radiation. In this field config- 
uration, therefore one has unknowingly imposed ingoing 
radiation to cancel the outgoing electromagnetic waves. 

Similarly, enforcing Kij — 7^ = might in part be 
thought of as implying the existence of ingoing gravi- 
tational radiation to cancel the outgoing gravity waves. 
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Nevertheless, in both cases, this remains a good approxi- 
mation to the physical system (with no ingoing wave) as 
long as the energy density contained in the radiation is 
small compared to the energy in orbital motion. 

Gravity waves enter in two ways: as estimated above 
there is an insignificant amount of "hidden" radiation in- 
duced by our choice of the CFC; there is also the emission 
of gravitational radiation by the orbiting binary system. 
The binary gravity- wave emission is estimated in our cal- 
culations by evaluating the multipole moments and using 
the appropriate formulas The fractional energy and 
angular momentum loss rate as determined by the multi- 
pole expansion method is quite small, e.g. J /ujJ ~ 10~ 4 
in all of our calculations Hence, it can be concluded 
that the energy in gravitational radiation is indeed small 
compared to the energy in orbital motion. 

The emission of gravity waves also induces a reaction 
force which we have incorporated into our hydrodynamic 
equations by the quadrupole formula. The radiation re- 
action force is so small, however, that it is difficult to 
discern it in the numerical results. In most of our cal- 
culations we simply neglect the back reaction terms and 
thereby obtain quasistatic orbit solutions. 



B. Solutions to Field Equations 

With a conformally flat metric, the constraint equa- 
tions for the field variables <f>, a, and (3 l reduce to simple 
Poisson-like equations in flat space. The Hamiltonian 
constraint equation |L7| for the conformal factor cf> be- 
comes 011, 



V 2 <b = -2nd> 5 



{l + U 2 )a-P+—K l3 K^ 



where a is the inertial mass-energy density 
a = p(l + e) + P , 



(10) 



(11) 



and p is the local proper baryon density which is simply 
related to the baryon number density n, p = pm^njN a-, 
where p is the mean molecular weight, the atomic 
mass unit, and Na is Avogadro's number, e denotes the 
internal energy per unit mass of the fluid, and P is the 
pressure. In analogy with special relativity we have also 
introduced a Lorentz-like variable 



1/2 



1 + WU, 



1/2 



1 • -' J i', I- j 



-1 1/2 



(12) 



where Ui is the spatial part of the covariant four veloc- 
ity. Here we explicitly write U 2 (in place of W 2 — 1 used 
in |l|-|3| ) because it emphasizes the extra velocity depen- 
dence here and in the equations of motion. 



In the Newtonian limit, the r.h.s. of Eq. ( |lC)| ) is dom- 
inated by the proper matter density p, but in rela- 
tivistic neutron stars there are also contributions from 
the internal energy density e, pressure P, and extrinsic 
curvature. This Poisson source is also enhanced by the 
generalized curved-space Lorentz factor (1 + U 2 ). This 
velocity factor becomes important as the orbit decays 
deeper into the gravitational potential and the orbital 
kinetic energy of the binary increases. 

It was pointed out in the appendix of || that in 
analogy to the velocity-dependent enhancement of the 
source for Eq. (|To|), the Poisson source for the v A post- 
Newtonian correction to the effective potential also ex- 
hibits velocity-dependence. This appendix has been mis- 
interpreted as a statement that we attribute the com- 
pression to a first post-Newtonian effect. We therefore 
wish to state clearly that the appendix in that paper 
was merely an illustration of how the effective grav- 
ity begins to show velocity dependence even in a post- 
Newtonian expansion. The velocity dependence of the 
post-Newtonian source is not the main compression driv- 
ing force. The compression derives mostly from the hy- 
drodynamic terms described herein. It is not obvious, 
however, at what post-Newtonian order the compression 
effect should be counted, since different authors have 
treated these terms differently. We return to this point 
below. 

In a similar manner the Hamiltonian constraint, 
together with the maximal slicing condition, provides an 
equation for the lapse function, 



V 2 (a</ 



-2p(l 



(3(C/ 2 + l)a 



3P 



16tt 



(13) 



Here again, the source is strengthened when the fluid is 
in motion through the presence of a U 2 + 1 factor and 



the KijK lJ term 



The momentum constraints JT^ ] provide an elliptic 
equation Q for the shift vector, 



J 3 > 



(14) 



p\= Uttf^-^'t^ + l)a 



1 dln{a/f) d d 2 d k 

dx> fe^ + fe^ _ S^fa^ ] ' (15) 

where we have introduced fl5f| the Lorentz contracted 
coordinate covariant momentum density, 



Si = oWU, 



(16) 



As noted previously and in Ref. Q , we only solve equa- 
tion ( |l4| ) for the small residual frame drag after the dom- 
inant 63 x r contribution to f3 has been subtracted. 
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III. RELATIVISTIC HYDRODYNAMICS 

The techniques of general relativistic hydrodynamics 
have been in place and well studied for over 25 years p5| . 
The basic physical processes which induce compression 
can be traced to completely general terms in the hydro- 
dynamic equations of motion. To illustrate this we first 
summarize the completely general derivation of the rel- 
ativistic covariant momentum equation in Eulerian form 
and identify the terms which we believe to be the domi- 
nant contributors to the relativistic compression effect. 

For hydrodynamic simulations it is convenient to ex- 
plicitly consider two different spatial velocity fields. One 
is Ui, the spatial components of the covariant four veloc- 
ity. The other is V* , the contravariant coordinate matter 
three velocity, which is related to the four velocity 



V 1 = 



IP ~ u l 



ft ■ 



(17) 



ft is convenient to select the shift vector ft such that the 
coordinate three velocity vanishes when averaged over 
the star, (V 1 ) — 0. This minimizes coordinate fluid mo- 
tion with respect to the shifting ADM grid. 
The perfect fluid energy-momentum tensor is 



V = oU^JJc + P 9i 



(18) 



However, it is convenient to derive the hydrodynamic 
equations of motion using the mixed form, 



(19) 



The vanishing of the spatial components of the diver- 
gence of the energy momentum tensor 



rp 







(20) 



leads to an evolution equation for the spatial components 
of the covariant four momentum, 



i d(s a ) 



°7 



Of 



1 d(SiVi-y) dP 
dx 1 



1 dcp? S a Sp 

2 dx 1 S l 



= 



(21) 



The covariant momentum equation is used because of 
its close similarity with Newtonian hydrodynamics. The 
first two terms are advection terms familiar from New- 
tonian fluid mechanics. The latter two terms are the 
pressure and gravitational forces, respectively. 

Expanding the gravitational acceleration into individ- 
ual terms we have 



7 1 d 

7 7 dxl 



adP 



S 



dft 



da 
dx % 



U 



, d In a 

dx % 



dx % 3 dx 
UjU k d^ k 
2 dx 1 



(22) 



Similar forms can be derived for the condition of 
baryon conservation and the evolution of internal energy 
0,0- However, the above momentum equation is suffi- 
cient for the present discussion. 

It is now worthwhile to consider the "gravitational" 
forces embedded in the expanded terms of Equation (|22]) . 
These result from the affine connection terms T^ X T^ in 
the covariant differentiation of T M ". 

The term containing da/dx l comes from the time-time 
component of the covariant derivative. It is of course the 
well known analog of the Newtonian gravitational force as 
can easily be seen in the Newtonian limit a — ► 1 — Gm/r. 

The term Sj(dft /dx 1 ) comes from the space-time co- 
variant derivative. In an orbiting system it is convenient 
to allow ft to follow the orbital motion of the stars. In 
our specific application || we let /? = w X R + P^lsid 
where u> is chosen to minimize matter motion on the 
grid. Hence, Co x R includes the major part of rota- 
tion plus frame drag. The quantity /3^™^ is the resid- 
ual frame drag after subtraction of rotation and is very 
small for the almost nonrotating stars which result from 
our calculations. With ft dominated by ui x R, the term 
S 3 (dft/dx 1 ) is predominantly a centrifugal force. 

The U dhxa/dx % term arises from the time-time com- 
ponent of the affine connection piece of the covariant 
derivative. The (UjUk/2)dft k /dx 1 term similarly arises 
from the space-space components. They do not have a 
Newtonian analog. As we shall see, these terms cancel 
when a frame can be chosen such that the whole fluid is 
at rest with respect to the observer (or in the flat space 
limit). However, for a star with fluid motion in curved 
space, they describe additional velocity-dependent forces. 

We identify the nonvanishing combination of these U 2 - 
dependent force terms and the Sj(dft/dx z ) term as the 
major contributors to the net compression driving force. 

This suggests some useful test problems for our hy- 
drodynamic simulations. For example, in simple uniform 
translation the effects of these terms must cancel to leave 
the stellar structure unchanged. Similarly, as discussed 
below, any fluid motion such that the four velocity can be 
taken as proportional to a simple Killing vector (e.g. rigid 
corotation) these force terms must cancel [ pips] ]. How- 
ever, for more general states of motion, e.g. noncorotating 
stars, differential rotation, meridional circulation, turbu- 
lent flow, etc., these forces do not obviously cancel, but 
must be evaluated numerically. 

Indeed, as discussed below, the sign of these terms is 
such that a lower energy configuration for the stars than 
that of rigid corotation can be obtained by allowing the 
fluid to respond to these forces. As we shall see, the 
numerical relaxation of binary stars from corotation (or 
any other initial spin configuration) produces a nonsyn- 
chronous (approximately irrotational) state of almost no 
intrinsic neutron-star spin in which the central density 
and gravitational binding energy increase. 
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A. Conformally Flat Relativistic Hydrodynamics 

The practical implementation of conformal flatness 
means that, given a distribution of mass and momentum 
on some manifold, we first solve the constraint equations 
of general relativity at each time for a given distribution 
of mass-energy. We then evolve the hydrodynamic equa- 
tions to the next time step. Thus, at each time slice we 
obtain a solution to the relativistic field equations and 
then can study the hydrodynamic response of the matter 
to these fields Q . 

For the CFC metric, the relativistic momentum equa- 



tion is derived by simply replacing -f jk 
Eq. M). 



dt 



■65, 



dim 



1 d 



dt 



dx l 



aaU' 



dlna 



<fS l V 3 + a 



dx % 



91n< 
) 

dx % 



dP 

dx % 

= 



s 



dp 
J 'dx 1 

(23) 



Here as in Eq. ([22|), the first term with da/dx 1 is the 
relativistic analog of the Newtonian gravitational force. 

In Eq. (^3|) there are two ways in which the effective 
gravitational force might increase for finite U 2 . One 
is that the matter contribution to the source densities 
for a or cf> are increased by factors of ~ 1 + U 2 [cf. 
Eqs. ([n]) and (p^)]. The more dominant effect, how- 
ever, is from the combination of the Sjd(3 J /dx l term and 
the U 2 [dlna/dx l - 2dln<p/dx i ] terms in Eq. (f|). 

As noted previously, these compression driving terms 
result from the affine connection part T^ a T mA of the co- 
variant differentiation of T^ v . These terms have no New- 
tonian analog but describe a general relativistic increase 
in the gravitational force as U 2 increases. As noted in 
|^|,D (see also Fig ^| below) for a binary, U 2 is approxi- 
mately uniform over the stars, and the increase in cen- 
tral density due to these additional forces scales as sw U 4 . 
This scaling, however, is the net result from a nontrivial 
cancellation of terms and must be treated carefully. We 
shall return to this point below. 

The proper way to determine the post-Newtonian or- 
der at which the compression driving terms enter would 
be to count the powers of c 2 which appear in the denom- 
inator of a term. For example, if we divide the last two 
terms in Eq. ( p3| ) by the gradient of the a term (the ana- 
log of the Newtonian gravitational force) we would obtain 
a ratio of order U 2 /c 2 which would be manifestly first 
post Newtonian. However, in the first post-Newtonian 
treatment of Wiseman j(J , these velocity terms were ex- 
plicitly disregarded. Thus, the effects of these terms 
could not have been present in that calculation. It is 
no surprise, therefore that no effect was observed in Ref. 

Also note that the 2d\n<f)/dx' 1 term in Eq. (|23j) en- 
ters with a sign such that the total C/ 2 -dependent con- 
tribution is further increased by about twice that from 
the dlna/dx 1 contribution alone. (The factor of 2 in 



front of the derivative comes from the requirement that 
4> 2 ~ (l/a) in the Newtonian limit [||.) 

A further increase of binding arises from the K 13 Kij 
terms in the field sources, but these terms are much 
smaller than the U 2 contributions for a binary system. 



B. Comment on the Relativistic Bernoulli Equation 

For comparison with other work in the literature it 
is instructive to discuss the derivation of the relativis- 
tic Bernoulli equation from equation (p2|). It has been 
pointed out (e.g. (l^]) that the hydrodynamics reduces 
to a simple equation for a fluid in which the velocity field 
can be represented by a Killing vector. In our notation 
this equation can be written, 



dm(f/' 



dP 



(24) 



The demonstration that the relativistic Bernoulli equa- 
tion (|J) is exactly reproduced from Eq. (|22|) and indeed 
for any case in which a Killing vector can be imposed, was 
recently brought to our attention by T. Nakamura [p9| . 
We summarize the derivation here in the conformally flat 
metric both for clarity and to show that conformal flat- 
ness does not violate this important constraint. 

To begin with, note that in the ADM formalism, the 
existence of a Killing vector is equivalent to being able to 
choose a shifted ADM grid such that V 1 = everywhere 
for the fluid. Next use Eq. ( fT^ ) to solve for f3 l and divide 
by a. The resulting equation for stationary motion is 



1 dP 

a dx % 



u, 



d 



dx i y^u* 



(aU 



L\2 



dlna 

dx l 



2U 



, <91n< 

dx % 



(25) 



The recovery of the relativistic Bernoulli equation re- 
quires that the r.h.s. = <91n[7* jdx 1 . With some straight- 
forward algebraic manipulation it is possible to show that 
all of the terms on the r.h.s. cancel except for one term 
from the /3 derivative, —<j)~ 4: U 2 dlnU t /dx' 1 . The comple- 
tion of the proof is simply to note that this term is equal 
to dlnUydx by Eq. (gjj). The result is Eq. (§J). 

It is instructive to consider the change in the relativis- 
tic Bernoulli equation when there is no Killing vector, 
i.e. Vj 0. Along the same lines of the derivation of 
Eq. ( p5[) , It can be shown [ffijf that the momentum equa- 
tion can be rewritten in our notation as, 



1 

aa 



7 1 d 
S t + S^ + - — (S l V 3 1 ) 
7 dx 3 



7 
1 dP 

a dx l 



IPUi 



d\n_ 

dx l 



dlnU 1 
dx l 



(26) 



The r.h.s. is just the relativistic Bernoulli equation in 
the limit that the l.h.s. vanishes. In general fluid flow, 
however, the l.h.s. contains not only the advection terms 
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(in brackets) , but also an additional surviving part of the 
P J derivative. 

It can be seen from this that imposing a Killing vector 
(V 1 — 0) means that only simple hydrostatic equilibrium 
is obtained for stationary systems. However, when non- 
trivial hydrodynamic motion is allowed, the extra forces 
embodied in the l.h.s. of Eq. ( p6| ) are manifest. The pres- 
ence of fluid motion not represented by a simple Killing 
vector, thus leads to a deviation from the simple rela- 
tivistic Bernoulli solution. Any attempt to model this 
deviation requires a careful treatment of the dynamical 
properties of the fluid described by the l.h.s. of Eq. (26). 



IV. CONSTRAINED HYDRODYNAMICS 

Further insight into the complexity of the physics con- 
tained in the relativistic equations of motion can be 
gained by considering some simple examples of con- 
strained hydrodynamics for which the answer is known. 
These pose useful tests of our numerical scheme. Since 
some have proposed that the effect we observe may be 
an artifact of numerical resolution or approximation, we 
present here a summary of various test problems designed 
to illustrate the stability of the numerics and also to 
compare with some of the calculations in the literature. 
These calculations also demonstrate that the compres- 
sion effect vanishes in the limiting cases which have been 
studied by others. Hence, they could not have been ob- 
served. They highlight the fact that the effect we ob- 
serve only appears in a strong field dynamic treatment 
which accounts for internal motion of stellar material in 
response to the binary and its effect on the star's self 
gravity. At present, ours may be the only existing result. 
This is consistent with the conclusion of |5(| based upon 
test particle dynamics. 



A. Bench Mark Calculations 

To test for the presence of the compression driving 
forces we consider two bench-mark initial calculations. 
The bench mark of no compression is that of an isolated 
star. In our three dimensional hydrodynamic calcula- 
tions, the single star structure is derived from Eq. J2^) 
in the limit 



Si = Ut = V 1 = (3 Z = 



(27) 



Hence, we again identify the dlna/dx 1 term with the 
relativistic analog of the Newtonian gravitational force. 
Eq. ( p8| ) also trivially reproduces to relativistic Bernoulli 
equation (f24|). 

We have of course tested our three dimensional cal- 
culations for single isolated stars. A single star remains 
stable on the grid indefinitely, except when baryon mass 
exceeds the maximum stable mass allowed by the TOV 
equations. Above the maximum TOV mass the stars be- 
gin to collapse on a dynamical timescale as they should. 
We have also checked that the grid resolution used in 
our binary calculations is adequate to produce the cor- 
rect central density, stellar radius, and gravitational mass 
of a single isolated star ||. Hence, it seems unlikely that 
inadequate grid resolution is the source of the compres- 
sion effect as some have proposed. 

In order to facilitate comparisons with the literature, 
and to avoid confusion over equation of state (EOS) is- 
sues, we have employed a simplistic r = 2 polytropic 
EOS, P = Kp r , where K = 1.8 X 10 5 erg cm 3 g~ 2 . 
This gives a maximum neutron-star mass of 1.82 Mq. 
The gravitational mass of a single tub = 1.625 Mq 
star in isolation is 1.51 Mq and the central density is 
p c = 5.84 x 10 14 g cm~ 3 . The compaction ratio is 
m/R — 0.15, similar to one of the stars considered in 
fll3f . Note that this EOS leads to stars with a lower com- 
paction ratio than the stars we considered in J2|,[| for 
which m/R as 0.2. Hence, the effects of tidal forces in 
the present calculations should be more evident. 

The bench mark in which compression is present is that 
of two equal mass stars in a binary computed with un- 
constrained hydrodynamics. The binary stars have the 
same baryon mass (tub = 1.625 Mq each), the same 
EOS, and a fixed angular momentum J = 2.5 x 10 11 cm 2 
(J/M B = 1.09 where Mb = 2ms). For these conditions 
the binary stars have U 2 — 0.025 and are at a coordinate 
separation of w 100 km. The stars are stable but close 
to the collapse instability. Hence, they have experienced 
some compression which has increased their central den- 
sity by 14% up to p c = 6.68 x 10 14 g cm" 3 . 

The central densities of these two bench marks are 
summarized in the first and last entries of Table Q. To 
compare with these bench-mark calculations we have 
computed equilibrium configurations for stars under the 
various conditions outlined below. The test for the pres- 
ence or absence of compression inducing forces will be the 
comparison of the numerically computed central density 
with that of a single isolated star or stars in a binary. 



The condition of hydrostatic equilibrium in isotropic co- 
ordinates is then trivially derived from Eq. (123) 



dP 

dx l 



d In a 

r ; 

dx 1 



(28) 



where the tilde denotes that the metric coefficients are 
evaluated in the fluid rest frame. The Newtonian limit 
of the right hand side is recovered as a — > 1 — Gm/r. 



B. Stars in Uniform Translation 

As a first nontrivial test, now consider a star as seen 
from an observer in an inertial frame which is in uniform 
translation with respect to the fluid. Choosing motion 
along the x coordinate, the fluid three velocity is, 
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JJX 

— = V x = Constant 



(29) 



However, the observer is still free to choose the ADM 
shift vector such that the computational grid remains 
centered on the star. That is, although Si,Ui,^ 0, we 
can still choose V 1 = 0. This gives a restriction on (3 l 



from Eq. (17), 



J xx aUx 
W 



(30) 



Note, that this is an ADM coordinate freedom. It is not 
equivalent to a coordinate boost. It is in fact a Killing 
vector which is convenient for numerical hydrodynamics. 
It allows the matter to remain centered on the grid even 
though the equations of motion are being solved for fluid 
which is not at rest with respect to the observer. 

With V % = 0, the x component of the momentum equa- 
tion (in equilibrium) becomes, 



dP 

dx 



S x d(3 x 
a dx 



(U 2 + 1) 



dlna U 3 U k dj jk 



dx 



With a CFC metric this becomes, 



dP_ _ S x d(3 x 
dx a dx 



(U 2 + 1) 



dh 



dx 



2U 2 



dx 



,d\nq 

dx 



(31) 



(32) 



There are now several differences between this expres- 
sion and that for an observer in the fluid rest frame. For 
one, there is the shift vector derivative d(3 x /dx. Even 
in uniform translation this derivative is nonzero due to 
the variations of the metric coefficients over the star (cf. 
Eq. ^). Also, the effective gravity is enhanced by the 
(U 2 + 1) velocity factor. The UjU k /2d^ jk /dx* term also 
appears. In addition, the effective source terms ( |Io| ) and 
( |l3| ) for the CFC metric coefficients are enhanced both 
by (U 2 + 1) factors and the K^K 1 ^ term. 

In spite of these differences, we nevertheless know that 
the locally determined pressure and inertial density must 
be the same as those determined for a star at rest. In- 
deed, since we can choose a Killing vector (V 1 = 0) these 
equations must reduce to the relativistic Bernoulli equa- 
tion d|). 

Thus, this is an important numerical test problem. 
We solve the full hydrodynamic equations explicitly, 
e.g. Eq. (p3[), under the initial condition of nonzero U x 
for a single star. The cancellations embedded in the hy- 
drodynamics are not obvious. Nevertheless, in the end, 
all of these effects must cancel to leave the stellar cen- 
tral density unchanged (except for a Lorentz contraction 
factor). 

To solve the uniform translation problem numerically 
we have applied the Hamiltonian and momentum con- 
straints to determine the metric coefficients. We then 
evolved the full hydrodynamic equations to equilibrium. 
Figure |^ shows the numerically evaluated central density 
for such translating stars as a function of U 2 . These stars 
were calculated with the EOS of M. This is compared 



with the central density for binary stars evolved at the 
same U 2 value using the same EOS. One can see that 
the translating stars maintain a constant central den- 
sity (within numerical error) as they should. In contrast, 
the central density of binary stars grows as « U A . This 
growth is the nontrivial net result from the velocity de- 
pendent terms in Eq. (p^). It is not obvious, however, to 
what post-Newtonian order this dependence corresponds. 

As summarized in Table [j], the central density for a 
uniformly translating r = 2 star with U 2 = 0.025 (for 
comparison with the binary bench-mark calculation) is 
5.90 x 10 14 g cm -3 . Within numerical accuracy, this 
central density is identical with that of an isolated star 
at rest. 

This is at least indicative that our observed growth in 
central density may not be numerical error as some have 
suggested (e.g. flC|Jl2] , |30|] ) . Such an error would likely 
be apparent in this test case. We argue that the differ- 
ence between simple translation and binary orbits relates 
to the physics of the binary system itself, in particular 
physics which is not apparent in uniform translation, an 
analysis of tidal forces, or a truncated expansion which 
does not contain sufficient terms to adequately describe 
the dynamical response of the fluid. 



C. Tidal Forces 



It has been pointed out jl 11 T^] that tidal forces are 
in the opposite sense to the compression driving forces 
discussed here. That is, tidal forces distort the stars 
and decrease the central density and therefore render the 
stars less susceptible to collapse. We have argued || 
that although such stabilizing forces are present in our 
calculations they are much smaller in magnitude than 
the velocity-dependent compression driving terms. Nev- 
ertheless, the evolution of the matter fields in a calcula- 
tion in which only tidal forces are present still represents 
a useful test of our numerical results. Stars in which 
only tidal forces act, should be stable and the central 
density should decrease rather than increase as the stars 
approach. 

To test the effects of tidal forces alone we have con- 
structed an artificial test calculation in which we place 
stars on the grid in a binary, but with no initial angular 
or linear momentum, i.e. J = and U 2 = 0. This ini- 
tial condition would normally evolve to an axisymmetric 
collision between the stars. However, after updating the 
matter fields, we artificially return the center of mass of 
the stars to the same fixed separation after each time 
step. We also reset to zero the mean velocity component 
directed along the line between centers. This sequence 
is repeated until the matter fields come to equilibrium. 
Since the velocity dependent forces eventually vanish, the 
only remaining forces are the pressure and static gravi- 
tational (including tidal) forces. 

Results as a function of separation distance are shown 
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in Table || for the T — 2 EOS and Table g for the re- 
alistic EOS used in [Q. For the realistic EOS the central 
density indeed decreases as the stars approach consistent 
with the expectations from Newtonian and relativistic 
tidal analyses For the T — 2 polytropic EOS, 

the central density also decreases as the stars approach 
and remain below the central density of an isolated star. 
The fact that this table is not monotonic at the inner- 
most point, however, may be due to a limitation of this 
numerical approximation for tidal forces as the ratio of 
separation to neutron-star radius diminishes. 

Although the tidal forces do indeed stabilize the stars, 
their effect on the central density is quite small (~ 0.2% 
decrease) compared to the net increase in density caused 
by the compression forces present for the binary. This is 
consistent with the relative order-of-magnitude estimates 
for these effects described in H). 



D. Stars in Rigid Corotation 

As a next nontrivial example, consider stars in a bi- 
nary system which are restricted to rigid corotation. In 
a recent series of papers, Baumgarte et al. have stud- 
ied neutron-star binaries using the same conformally flat 
metric. Their work differs from ours in that rather than 
solving the hydrodynamic equations, they describe the 
four velocity field by a Killing vector whereby the stars 
are forced to corotate rigidly. They also impose spatial 
symmetry in the three Cartesian coordinate planes so 
that they can solve the problem in only one octant. One 
should keep in mind, however, that rigid corotation is not 
necessarily the lowest energy configuration or the most 
natural |32f | final state for two neutron stars near their 
final orbits. This assumption, though artificial, is never- 
theless a means to constrain and simplify the fluid motion 
degrees of freedom. It is much easier to implement and 
therefore becomes an interesting test problem for codes 
seeking to explore the true hydrodynamic evolution of 
close binaries. 

Indeed, it is possible to show [^U that in this limit, 
the neutron star hydrostatic equilibrium can be described 
by a simple Bernoulli equation in which the compression 
driving force terms are absent except for a weak velocity 
dependence. Analytically, the reason for this is trivially 
obvious from Eq. (^6|). The existence of a Killing vector 
is equivalent to setting V % = globally. Choosing the 
ADM coordinates to remain centered on the stars, in 
steady state the time derivatives vanish along with the 
rest of the l.h.s of Eq. (p6|). Only the relativistic Bernoulli 
equation ( p4[ ) survives. 

It is not surprising, therefore that in |D| it has been 
demonstrated that in this special symmetry, the central 
density of the stars does not increase (within numerical 
error) as the stars approach the inner most stable circular 
orbit. In very close orbits the density actually decreases 
relative to the central density of stars at large separation. 



They also find that the orbit frequency remains close to 
the Newtonian frequency. Both of these results are inter- 
esting in that they confirm that the compression effect 
does not occur (as it should not) in this special symme- 
try. They also demonstrate that conformal flatness is not 
the source of the compression. 

Accepting the results of |l3| as correct, this then be- 
comes another important test of our calculations. That 
is, if we artificially impose rigid corotation, then the cen- 
tral density should remain nearly constant until the stars 
are close enough that tidal effects cause the central den- 
sity to decrease rather than increase. 

Imposing rigid corotation, however, is not a trivial 
test problem to implement without completely replac- 
ing the hydrodynamic equations with the corresponding 
Bernoulli solution of (Indeed, we have done this 

ifjlf and reproduce the results of Jl3| quite well). More- 
over, we have found that directly modifying the hydrody- 
namic equations in an attempt to mimic a dynamically 
unstable configuration is difficult. One might think that 
the simplest way to implement corotation would be to im- 
pose a high fluid viscosity. Indeed high viscosity would 
resist the hydrodynamic forces described herein. How- 
ever, a high fluid viscosity also resists the much weaker 
tidal forces and prevents the numerical relaxation to qua- 
sistatic equilibrium. It is thus difficult to achieve tidal 
locking by simply increasing the viscosity. 

Instead, we introduce artificial forces on the fluid which 
continually drive the system toward a state of rigid coro- 
tation while allowing the system to at least somewhat 
respond hydrodynamically. To do this we define acceler- 
ations (t/^FUgid necessary to achieve rigid rotation by 



(^i)Rigid = ^7 



(33) 



where f7, are components of the rigidly corotating co- 
variant four velocity in the x — y orbit plane. These are 
determined by requiring that (3 l — (u> x r) 1 and setting 
V i = in Eq. @. 



U„ = 



U x 



-uycfr 



(34) 



(35) 



where R is the coordinate distance from the center of 
mass of the binary. 

At each time step we then update the momentum den- 
sity using a combination of the hydrodynamic and coro- 
tating acceleration terms, 



Ui = /(C/i)Rigid + (1 - /XCOHydr 



(36) 



where (Ui)nydro is the acceleration from the full hydro- 
dynamic equation of motion [Eq. (p3|)]. 
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Numerically, we find that if / is small (< 0.2) the hy- 
drodynamic forces dominate and corotation is not ob- 
tained. On the other hand, for / > 0.2 the system is 
not stable, i.e. the stars deform and the velocities be- 
come erratic. We have therefore run with / = 0.2 which 
temporarily produces a velocity field which is close to 
rigid corotation. That is, the residual three velocities are 
damped to a fraction of the orbit speed. This is, perhaps, 
good enough to make qualitative comparisons with the 
expectations from a truly corotating system. 

Starting from the unconstrained initial configuration, 
we find that when the stars have achieved approximate 
corotation the central density has decreased from 6.68 x 
10 14 g cm -3 to 5.90 x 10 14 g cm 3 which is close to 
the value for stars in isolation (5.84 x 10 14 g cm -3 ). The 
calculated gravitational mass is slightly greater than that 
of the unconstrained binary. However, with the large 
artificial force terms needed to approximate corotation, 
gravitational mass is not a well defined quantity in this 
simulation. Also, the orbit frequency was not sufficiently 
converged for a meaningful comparison. 

E. The Spin of Binary Stars 

As noted above our simulations indicate that neutron 
stars relax to a state of almost no intrinsic spin. In a 
separate paper [[33| we analyze the nature and formation 
of this state in more detail. For the present discussion, 
however, we summarize in Figure ^ a study of the relax- 
ation to this state from states of arbitrary initial rigid 
rotation (including corotation). 

As a means to distinguish the intrinsic spin motion of 
the fluid with respect to a non-orbiting distant observer, 
we define a quantity which is analogous to volume aver- 
aged intrinsic stellar spin in the orbit plane, 




where yi, Zi = 0) is the coordinate center of mass of 
each star. 

In this study we have imposed an initial angular ve- 
locity ujs in the corotating frame to obtain various initial 
rigidly rotating spin angular momenta (including coro- 
tation, uj s = 0), but for fixed total J/M% = 1.4. We 
have considered spin angular frequencies in the range 
—900 < ujs < 900 rad sec -1 , corresponding to —0.03 < 
Js/m% < 0.17. We then let the system evolve hydrody- 
namically with the stars maintained at zero temperature. 

In Ref. @ we showed that neutrino emission is suffi- 
cient to radiate away the released gravitational energy 
and keep the stars at near zero temperature until just 
before collapse. This is the reason that we have treated 
this as a relaxation problem. That is, unlike a true hy- 
drodynamic calculation, the relaxation calculation pre- 
sented here, assumes that the stars radiate efficiently and 



stay at zero temperature. Therefore, this evolution does 
not need to conserve energy or circulation. This relax- 
ation assumption is the reason the stars can evolve to a 
different spin (lower energy) state without violating the 
circulation theorem. 

Figure [| shows the spin Js /m 2 B as a function of time 
for each initial condition. In each case, the system re- 
laxed to a state of almost no net spin within about three 
sound crossing times (t ~ 0.6 msec). These calculations 
suggest that rapidly spinning neutron stars in close orbits 
are unstable. The true evolution time, however, would 
be much longer. 

We also note that the quantity J a[y/l + U 2 — l]dV 
decreased as the system evolved from rigid rotation to 
hydrodynamic equilibrium. Since this quantity is related 
to the kinetic energy of the binary, this indicates that the 
hydrodynamic lowest energy state is one of lower kinetic 
energy (for fixed total angular momentum) than that of 
rigid rotation. 

As far as the compression effect is concerned, one 
wishes to know whether the response of the stars is sim- 
ply due to that fact that they have no spin (and therefore 
no internal centrifugal force to support them against the 
compression forces), or whether more complex fluid mo- 
tion within the star itself affects the stability. To test 
this, we have constructed stars of no spin (Jg = 0) by 
simply damping the residual motion to that of Js = 
after each update of the velocity fields. 

Since this no-spin state is so close to the true hydro- 
dynamic equilibrium, this produced stable Js = equi- 
librium stars for the binary. For this case, the central 
density converges to p c = 6.56 x 10 14 g cm -3 which is 
very close to the high value for the unconstrained hydro- 
dynamics. This result would seem to indicate that most 
of the increase in density can be attributed to the veloc- 
ity with respect to the corotating frame generated by the 
fact that the stars have almost no spin. 

V. DISCUSSION 

For clarity, we summarize in this section our conclu- 
sions regarding why the neutron-star compression effect 
was not observed in some other recent works. 

First consider post-Newtonian expansions. In the work 
of Wiseman Q the force terms containing U 2 were explic- 
itly deleted from the computation of the stellar structure 
[cf. Eq. (8) in that paper]. Only the dlna/dx term was 
included. The recovery of simple hydrostatic equilibrium 
was thus unavoidable. 

The PN orbiting ellipsoids of Shibata et al. included 
more terms. Indeed, it was noted that there are two ef- 
fects at 1PN order. One is the self gravity of each star 
of the binary and the other is the gravity acting between 
the stars. In their calculations the self gravity dominates 
causing the stars to become more compact. This is con- 
sistent with the compression effect described here in the 
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sense that relativistic corrections can dominate over New- 
tonian tidal forces. However, the self gravity terms in J?]] 
appear to only include the usual 1PN terms which would 
equally apply to stars in isolation. Hence, the velocity- 
dependent compression driving terms are probably not 
present. 

Their results for stars in corotation are consistent with 
ours under the same constraint. They also note that ap- 
proaches in which PN corrections to the gravity between 
the stars is included without also including the correc- 
tions to the self gravity (as in |9|) can be misleading. 

In the work of Lombardi et al. || both corotating and 
irrotational equilibria were computed. However, in their 
calculations it appears that the stars become less com- 
pact as they approach contrary to our results and the 
results of M. It may be that the reason for this is that 
in Lombardi et al. the post-Newtonian corrections to self 
gravity were only computed for stars " instantaneously at 
rest" . The authors chose to " exclude the spin kinetic en- 
ergy contribution to the self energy". It is such terms, 
however, which we identify with the compression effect. 

The conformally flat corotating equilibria computed 
by Baumgarte et al. jl3| are consistent with our results. 
Since their stars were restricted to rigid corotation, only 
the hydrostatic Bernoulli solution would result. They 
could not have observed the compression forces which 
result from fluid motion with respect to the corotating 
frame. 

We have argued in this paper that if one wishes to ex- 
plore this effect, it would be best to apply a complete un- 
constrained strong-field relativistic hydrodynamic treat- 
ment for stars which are not in corotation. In this re- 
gard, a recent paper p4[ has come to our attention in 
which hydrodynamic simulations of both corotating and 
irrotational binaries have been studied in a first post- 
newtonian approximation to conformally-flat gravity but 
using the full relativistic hydrodynamics equation (22). 
For both corotational and irrotational stars the central 
density is observed to oscillate about a value which is less 
than that of isolated stars. Hence, the authors conclude 
that no compression effect is present. 

Since this calculation contains many of the higher or- 
der terms to which we attribute the compression effect, 
it is not immediately obvious why the compression effect 
was not observed. This may indeed be a real contradic- 
tion. We suggest, however, that this simulation did not 
observe the effect because of their use of an unrealistically 
soft r = 1.4 EOS. The authors chose this EOS because 
the stars become so extended that one can compute arbi- 
trarily close binaries without encountering the relativistic 
inner orbit instability. For the irrotational stars (model 
Be in [ll|l ) , which is the only simulation that might have 
observed the compression effect, the compaction ratio is 
only M/R = 0.023. Hence, a 1.45 M Q neutron star would 
have an unrealistic radius of 93 km. 

However, since they have simulated very extended 
stars at very close separation, the tidal forces are much 
stronger relative to the relativistic compression driving 



terms than in any of the simulations which we have done. 

The ratio of the stabilizing tidal correction AE t idai 
to the destabilizing energy from compression AE cornp 
should scale 0,01 as 



AE K I 7 



(38) 



where R is the neutron star radius and r is the orbital 
separation. For model Be in Jlij we estimate that this 
ratio is ^200 times greater than any of the binary stars 
we have considered. Hence, it is quite likely that the au- 
thors have simply chosen an unrealisticly soft equation of 
state for which the tidal forces dominate over compres- 
sion. It might be very interesting to see the results from 
a similar study for stars with a realistic compaction ratio 
and several radii apart. 

Concerning tidal expansions, in Brady & Hughes 0] 
an attempt was made to analyze the stability of a central 
star perturbed by an orbiting point particle. The met- 
ric and stress-energy were perturbed in terms of order 
e = fi/R where /i is the point particle mass and R its 
coordinate distance from the central star. The Einstein 
equation was then linearized to terms of order e. The 
result of this linearization was that the only possible cor- 
rection to the central density was a single monopole term 
of order [ijR ~ v 2 . However, in our numerical results as 
shown in Figure |^. the central density is observed to in- 
crease as v A . Hence, it may be that the expansion of Ref. 
]iot was truncated at too low order to observe the com- 
pression effect described here. The main reason that they 
could not observe the effect, however, is that the terms 
involving motion of the central star were discarded. We 
attribute the compression effect to an enhancement of the 
self gravity due to motion of the stars with respect to the 
corotating frame. Hence, the neglect of terms involving 
motion of the central star precludes the possibility of ob- 
serving the effect. 

We believe that the same conclusion is true in the 
treatments by Refs. |ll]|l^]. The analysis of Flanagan 
is based upon the method of matched asymptotic 
expansion. The metric is approximated 



9^ = Vitv + + h 



(39) 



where the superscript NS refers to the self contribution 
from one star and B refers to the contribution from a 
distant companion. The internal gravity of a static neu- 
tron star hp„ is expanded to all orders. The binary tidal 
contribution h^ v is expanded in powers of the ratio of 
stellar radius to orbital separation. 

First we suggest that such a decomposition may be 
questionable for a close neutron-star binary. In our met- 
ric one can write the metric perturbation as 



h l3 = (0 4 - 1)^ 



(40) 



The conformal factor </> is a solution to a Poisson equation 
involving source terms from the two stars. Between the 



11 



stars, the only source of the fields arises from the K^K 1 ^ 
terms which are quite small. Hence, neglecting KijK^ 
terms, (j) is additive in the "vacuum" between the stars, 



= 1 



m 2 



2\r-r x \ 2\r-r 2 \ 



(41) 



Expanding hij around star 1 in the presence of a distant 
companion 2 we have 



mi 



m 2 



r-i | \r-r 2 \ 



6 



m 2 



4 \|r- ri | \r - r 2 \ 



= h 



NS 

fiV 



h.„. + cross terms 



(42) 



However, for the binary systems we have considered, the 
cross terms are ~ 15% to 20% of the sum + h^ v . 
Hence, they can not be neglected. The errors associated 
with this decomposition may be part of the reason that 
the compression effects are not apparent in this work. 

A related concern is with the expansion of the stress- 
energy tensor in jll]]. We have noted that most of the 
compression arises from the net effect of velocity depen- 
dent terms in the covariant derivative of the stress-energy 
tensor. In [|ll| the stress energy is expanded is pow- 
ers of the curvature R~ m . The author states Q "We 

(2) — 

assume initial conditions of vanishing Tp, v , so that the 
only source for perturbations is the external tidal field." 
An analysis which only considers perturbations from the 
external tidal field (and not motions of the fluid) will 
not observe the compression effect. The result of Jll| is 
that the central density is unchanged until tidal forces 
enter at 0(R 6 ). This is consistent with our results in the 
limit of only tidal perturbations acting on the stars. It 
is not clearto us, however, to what degree the velocity 
dependent terms are included or excluded by this expan- 
sion. A more careful recent revision (E. Flanagan, Priv. 
Comm.) shows an effect coming in a lower order, but not 
necessarily as strong as we have noted. 

In the paper of Thorne jl2| , a similar tidal expansion 
is applied. In that work only the stabilizing effect of 
tidal forces was considered along with the stabilizing ef- 
fect of rotation. However, the increased self gravity from 
velocity-dependent forces was not included. Hence, the 
conclusions of |1^| are consistent with our results based 
upon tidal forces. So are the Newtonian tidal effects com- 
puted in |] . 



VI. CONCLUSIONS 

The results of this study (cf. Table |) are that we see 
almost no difference between the central density of an 
isolated star and a binary star in which rigid corotation 
has been artificially imposed, or one in which only tidal 
effects are included. Indeed, in the case of tidal forces 



alone, the central density in our simulations actually de- 
creases as stars approach, consistent with other works. 

An increase in the central density is only apparent in 
our binary simulations for stars with fluid motion with 
respect to the corotating frame. (Specifically we consid- 
ered stars of low intrinsic spin in a binary.) In such cases 
there is no simple Killing vector which can be imposed to 
cancel the compression driving forces. We have argued 
here and in || that the main compression effect arises 
from the net result of velocity-dependent hydrodynamic 
terms p4| . These terms arise from the affine connection 
part of the covariant differentiation of the stress-energy 
tensor. 

We show here that the compression effect would not 
have been observed in a study of tidal forces or any model 
which artificially imposes rigid corotation of the fluid. A 
proper treatment must consider all of the force terms ap- 
parent in the momentum equation ( j2^ ) to sufficient order 
that their effects on the fluid self gravity survive. A simi- 
lar conclusion has been reached in |3(| based on test par- 
ticle dynamics near a Schwarzschild black hole. In that 
work it is concluded that at least 2.5 post-Newtonian par- 
ticle dynamics is necessary before a dynamical collapse 
instability is manifest. 

We argue that the results of this study are thus con- 
sistent with results in a number of recent papers |4j-|l5|l 
which have analyzed the stability of binary stars in var- 
ious approximations and limits and see no effect. Since 
we do not disagree with the lack of a compression effect 
in the limits which they have imposed, we conclude that 
the existence or absence of the neutron-star compression 
effect has not yet been independently tested. 

Therefore, if one wishes to explore this effect, it would 
be best to apply a complete unconstrained strong-field 
relativistic hydrodynamic treatment employing an EOS 
which produces realistically compact neutron stars. An- 
other alternative, however, might be to study the quasi- 
equilibrium structure of nonspinning irrotational binary 
stars at sufficiently high order. In this regard a recently 
proposed formalism [p5[ to compute quasi-equilibria for 
nonsynchronous binaries may be of some use. We have 
begun calculations in this independent formalism. The 
results will be reported in a forthcoming paper. 

Regarding the existence of this low spin state, we find 
that such a state represents the unconstrained hydrody- 
namic equilibrium for a close binary. In Newtonian the- 
ory, stars are driven to corotation by tidal forces. How- 
ever in H| it has been shown that Newtonian tidal forces 
are insufficient to produce corotation before neutron-star 
merger unless the viscosity is unrealistically high. Nev- 
ertheless, in the absence of strong tidal forces, neutron 
stars stars gradually spin down. Therefore, even apart 
from the hydrodynamic effects described here, stars of 
low spin are likely to be members of close binaries. The 
hydrodynamic effects described herein, however, could 
hasten the spin down as stars approach their final orbits 
and cause the stars to become more compact. 
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FIG. 1. The scaled Cotton- York tensor component (j" r m~ 
as a function of proper radius r/m for a maximally rotating 
a = m Kerr black hole. This quantity is a measure of the 
deviation from conformal flatness. 



FIG. 2. Numerically evaluated central density for a uni- 
formly translating star (lower curve) as a function of 
U 2 = U l Ui. This is compared with the central density for 
binary stars (upper curve) with the same average U 2 value. 
These calculations utilized the EOS of Ref. [2]. 



FIG. 3. Intrinsic neutron star spin Js/m 2 B as a function 
of coordinate time. The curves are labeled by the initial an- 
gular velocity los (in units of 100 rad sec - ) relative to the 
corotating frame. 
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TABLE I. Central density for m s 
various conditions using a T = 2 EOS. 



Environment 



Constraints 



1.625 Mq stars in 



p c (10 14 g cm-*) 



Single Star 
Single Star 
Binary 
Binary 
Binary 
Binary 



Hydrostatic 
Uniform Translation 
Tidal Only 
Rigid Corotation 
Rigid No Spin 
Full Hydrodynamics 



5.84 
5.90 
5.82 
5.90 
6.56 
6.68 



TABLE II. Central density vs. coordinate separation be- 
tween centers for = 1.625 M Q (T = 2) stars in which only 
tidal forces are included. The neutron star radius (in isotropic 
coordinates) is 12 km. 



Separation (km) 



p c (10 14 g cm" 3 ) 



41.8 
50.6 
81.8 

oo 



5.821 
5.816 
5.821 
5.837 



TABLE III. Same as Table II but for the EOS of ref. [2] 
The neutron star radius (in isotropic coordinates) is 6 km. 



Separation (km) 



p c (10 14 g cm" 3 ) 



31.2 
37.4 
64.8 
103.8 

oo 



14.15 
14.16 
14.20 
14.25 
14.30 
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